Root_Atlas.rds (or get it by running through notebook 2, 3, 4, 5, 6 & 7)
Root_Atlas_spliced_unspliced_raw_counts.rds
rm(list=ls())
# Set the working directory to where folders named after the samples are located.
# The folder contains spliced.mtx, unspliced.mtx, barcodes and gene id files, and json files produced by scKB that documents the sequencing stats.
setwd("/scratch/AG_Ohler/CheWei/scKB")
# Load libraries
suppressMessages(library(Seurat))
suppressMessages(library(CytoTRACE))
suppressMessages(library(RColorBrewer))
suppressMessages(library(ggplot2))
suppressMessages(library(grid))
Warning message: “The ScanoramaCT python module is not accessible. The iCytoTRACE function for integration across multiple datasets will be disabled. Please follow the instructions in https://github.com/gunsagargulati/CytoTRACE to install the necessary Python packages for this application.”
sessionInfo()
R version 4.0.3 (2020-10-10) Platform: x86_64-conda-linux-gnu (64-bit) Running under: CentOS Linux 7 (Core) Matrix products: default BLAS/LAPACK: /fast/home/c/chsu/anaconda3/envs/seu3/lib/libopenblasp-r0.3.12.so locale: [1] LC_CTYPE=en_US.utf-8 LC_NUMERIC=C [3] LC_TIME=en_US.utf-8 LC_COLLATE=en_US.utf-8 [5] LC_MONETARY=en_US.utf-8 LC_MESSAGES=en_US.utf-8 [7] LC_PAPER=en_US.utf-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.utf-8 LC_IDENTIFICATION=C attached base packages: [1] grid stats graphics grDevices utils datasets methods [8] base other attached packages: [1] ggplot2_3.3.2 RColorBrewer_1.1-2 CytoTRACE_0.1.0 Seurat_3.1.5 loaded via a namespace (and not attached): [1] httr_1.4.2 tidyr_1.1.2 jsonlite_1.7.1 [4] viridisLite_0.3.0 splines_4.0.3 leiden_0.3.5 [7] ggrepel_0.8.2 globals_0.13.1 pillar_1.4.6 [10] lattice_0.20-41 glue_1.4.2 reticulate_1.18 [13] uuid_0.1-4 digest_0.6.27 colorspace_2.0-0 [16] cowplot_1.1.0 htmltools_0.5.0 Matrix_1.2-18 [19] plyr_1.8.6 pkgconfig_2.0.3 tsne_0.1-3 [22] listenv_0.8.0 purrr_0.3.4 patchwork_1.1.0 [25] scales_1.1.1 RANN_2.6.1 Rtsne_0.15 [28] tibble_3.0.4 generics_0.1.0 ellipsis_0.3.1 [31] withr_2.3.0 repr_1.1.0 ROCR_1.0-11 [34] pbapply_1.4-3 lazyeval_0.2.2 survival_3.2-7 [37] magrittr_1.5 crayon_1.3.4 evaluate_0.14 [40] future_1.20.1 parallelly_1.21.0 nlme_3.1-150 [43] MASS_7.3-53 ica_1.0-2 tools_4.0.3 [46] fitdistrplus_1.1-1 data.table_1.13.2 matrixStats_0.57.0 [49] lifecycle_0.2.0 stringr_1.4.0 plotly_4.9.2.1 [52] munsell_0.5.0 cluster_2.1.0 irlba_2.3.3 [55] compiler_4.0.3 rsvd_1.0.3 rlang_0.4.8 [58] ggridges_0.5.2 pbdZMQ_0.3-3.1 IRkernel_1.1.1.9000 [61] rappdirs_0.3.1 RcppAnnoy_0.0.16 htmlwidgets_1.5.2 [64] igraph_1.2.6 base64enc_0.1-3 gtable_0.3.0 [67] codetools_0.2-18 reshape2_1.4.4 R6_2.5.0 [70] gridExtra_2.3 zoo_1.8-8 dplyr_1.0.2 [73] uwot_0.1.8 future.apply_1.6.0 KernSmooth_2.23-18 [76] ape_5.4-1 stringi_1.5.3 parallel_4.0.3 [79] IRdisplay_0.7.0 Rcpp_1.0.5 sctransform_0.3.1 [82] vctrs_0.3.4 png_0.1-7 tidyselect_1.1.0 [85] lmtest_0.9-38
# Plotting function for cell tyepes and time zone
plot_anno <- function(rc.integrated){
order <- c("Quiescent Center", "Ground Tissue","Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Metaphloem & Companion Cell","Protophloem", "Xylem", "Procambium","Xylem Pole Pericycle","Phloem Pole Pericycle", "Protoxylem", "Metaxylem", "Unknown")
palette <- c("#9400d3", "#DCD0FF","#5ab953", "#bfef45", "#008080", "#21B6A8", "#82b6ff", "#0000FF","#e6194b", "#dd77ec", "#9a6324", "#ffe119", "#ff9900", "#ffd4e3", "#9a6324", "#ddaa6f", "#EEEEEE")
rc.integrated$celltype.anno <- factor(rc.integrated$celltype.anno, levels = order[sort(match(unique(rc.integrated$celltype.anno),order))])
color <- palette[sort(match(unique(rc.integrated$celltype.anno),order))]
p1 <- DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno", cols=color)
p2 <- DimPlot(rc.integrated, reduction = "umap", group.by = "time.anno",order = c("Distal Columella","Proximal Columella","Distal Lateral Root Cap","Proximal Lateral Root Cap","Maturation","Elongation","Meristem"),cols = c("#e5f5e0", "#a1d99b", "#31a354", "#deebf7", "#3182bd", "#fee0d2","#de2d26"))
options(repr.plot.width=20, repr.plot.height=8)
gl <- lapply(list(p1, p2), ggplotGrob)
gwidth <- do.call(unit.pmax, lapply(gl, "[[", "widths"))
gl <- lapply(gl, "[[<-", "widths", value = gwidth)
gridExtra::grid.arrange(grobs=gl, ncol=2)
}
zscore <- function(x){(x-mean(x))/sd(x)}
# Read in atlas
rc.integrated <- readRDS('./Root_Atlas.rds')
rc.su.counts <- readRDS('./Root_Atlas_spliced_unspliced_raw_counts.rds')
table(rc.integrated$celltype.anno)
Quiescent Center Columella
158 8535
Lateral Root Cap Atrichoblast
18396 13380
Trichoblast Cortex
12361 11073
Endodermis Metaphloem & Companion Cell
11369 3224
Protophloem Procambium
939 9307
Xylem Pole Pericycle Phloem Pole Pericycle
12213 4634
Protoxylem Metaxylem
3195 1643
# Index for ground tissue trajectory
end.cor.traj.idx <- which(rc.integrated$celltype.anno == "Endodermis" | rc.integrated$celltype.anno == "Cortex" | rc.integrated$celltype.anno == "Quiescent Center")
# Extract ground tissue
end.cor.integrated <- subset(rc.integrated, cells = colnames(rc.integrated)[end.cor.traj.idx])
# Run UMAP
end.cor.integrated <- RunUMAP(end.cor.integrated, reduction = "pca", dims = 1:8, umap.method = "umap-learn", metric = "correlation")
end.cor.integrated[["umap"]]@cell.embeddings[,1] <- end.cor.integrated[["umap"]]@cell.embeddings[,1]*-1
end.cor.integrated[["umap"]]@cell.embeddings[,2] <- end.cor.integrated[["umap"]]@cell.embeddings[,2]*-1
#u2 <- rc.integrated@reductions$umap@cell.embeddings[,1]
#u1 <- rc.integrated@reductions$umap@cell.embeddings[,2]
#rc.integrated@reductions$umap@cell.embeddings[,1] <- u1
#rc.integrated@reductions$umap@cell.embeddings[,2] <- u2
plot_anno(end.cor.integrated)
# Prepare expression matrix for CytoTRACE
expression_matrix <- end.cor.integrated@assays$integrated@data
expression_matrix[which(expression_matrix < 0)]=0
expression_matrix <- as(expression_matrix, "dgCMatrix")
end.cor.integrated@assays$integrated@counts <- expression_matrix
# Run CytoTRACE
results <- CytoTRACE(as.matrix(end.cor.integrated@assays$integrated@counts), ncores = 16, subsamplesize = 1000)
end.cor.integrated$CytoTRACE <- results$CytoTRACE
The number of cells in your dataset exceeds 3,000. CytoTRACE will now be run in fast mode (see documentation). You can multi-thread this run using the 'ncores' flag. To disable fast mode, please indicate 'enableFast = FALSE'. CytoTRACE will be run on 23 sub-sample(s) of approximately 983 cells each using 16 / 16 core(s) Pre-processing data and generating similarity matrix... Calculating gene counts signature... Smoothing values with NNLS regression and diffusion... Calculating genes associated with CytoTRACE... Done
options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(end.cor.integrated, features = "CytoTRACE", pt.size=0.5)+ scale_colour_gradientn(colours = rev(brewer.pal(11,"Spectral")))
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
# Save Seurat object
saveRDS(end.cor.integrated,'./supp_data/Ground_Tissue_Atlas.rds')
# Prepare scVelo input
rc.su.counts <- subset(rc.su.counts, cells = colnames(rc.su.counts)[end.cor.traj.idx])
sr <- rc.su.counts@assays$spliced_RNA@counts
sr <- sr[rownames(end.cor.integrated@assays$integrated@data),]
ur <- rc.su.counts@assays$unspliced_RNA@counts
ur <- ur[rownames(end.cor.integrated@assays$integrated@data),]
ar <- end.cor.integrated@assays$RNA@counts
ar <- ar[rownames(end.cor.integrated@assays$integrated@data),]
sr <- sr/ar;
ur <- ur/ar;
sr <- as.matrix(sr)
ur <- as.matrix(ur)
sr[is.nan(sr)] = 0;
ur[is.nan(ur)] = 0;
colnames(sr) <- colnames(end.cor.integrated)
colnames(ur) <- colnames(end.cor.integrated)
int <- as.matrix(end.cor.integrated@assays$integrated@counts)
spliced <- int*sr;
unspliced <- int*ur;
sg <- intersect(rownames(spliced), rownames(unspliced));
spliced <- spliced[match(sg, rownames(spliced)),];
unspliced <- unspliced[match(sg, rownames(unspliced)),];
meta <- end.cor.integrated@meta.data[,grep("time.anno|celltype.anno|time.celltype.anno|CytoTRACE",colnames(end.cor.integrated@meta.data))];
var <- sg;
pca_int <- end.cor.integrated@reductions$pca@cell.embeddings;
umap_int <- end.cor.integrated@reductions$umap@cell.embeddings;
save(spliced, unspliced, meta, var, pca_int, umap_int, file = "./supp_data/Ground_Tissue_Atlas_scVelo_input.RData")
# Index for epidermis + LRC trajectory
end.cor.traj.idx <- which(rc.integrated$celltype.anno == "Quiescent Center" | rc.integrated$celltype.anno == "Atrichoblast" | rc.integrated$celltype.anno == "Trichoblast" | rc.integrated$celltype.anno == "Lateral Root Cap" );
# Extract epidermis + LRC
end.cor.integrated <- subset(rc.integrated, cells = colnames(rc.integrated)[end.cor.traj.idx])
# Run UMAP
end.cor.integrated <- RunUMAP(end.cor.integrated, reduction = "pca", dims = 1:8, umap.method = "umap-learn", metric = "correlation")
plot_anno(end.cor.integrated)
# Prepare expression matrix for CytoTRACE
expression_matrix <- end.cor.integrated@assays$integrated@data
expression_matrix[which(expression_matrix < 0)]=0
expression_matrix <- as(expression_matrix, "dgCMatrix")
end.cor.integrated@assays$integrated@counts <- expression_matrix
# Run CytoTRACE
results <- CytoTRACE(as.matrix(end.cor.integrated@assays$integrated@counts), ncores = 16, subsamplesize = 1000)
end.cor.integrated$CytoTRACE <- results$CytoTRACE
The number of cells in your dataset exceeds 3,000. CytoTRACE will now be run in fast mode (see documentation). You can multi-thread this run using the 'ncores' flag. To disable fast mode, please indicate 'enableFast = FALSE'. CytoTRACE will be run on 44 sub-sample(s) of approximately 1007 cells each using 16 / 16 core(s) Pre-processing data and generating similarity matrix... Calculating gene counts signature... Smoothing values with NNLS regression and diffusion... Calculating genes associated with CytoTRACE... Done
options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(end.cor.integrated, features = "CytoTRACE", pt.size=0.5)+ scale_colour_gradientn(colours = rev(brewer.pal(11,"Spectral")))
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
saveRDS(end.cor.integrated,'./supp_data/Epidermis_LRC_Atlas.rds')
# Prepare scVelo input
rc.su.counts <- readRDS('./Root_Atlas_spliced_unspliced_raw_counts.rds')
rc.su.counts <- subset(rc.su.counts, cells = colnames(rc.su.counts)[end.cor.traj.idx])
sr <- rc.su.counts@assays$spliced_RNA@counts
sr <- sr[rownames(end.cor.integrated@assays$integrated@data),]
ur <- rc.su.counts@assays$unspliced_RNA@counts
ur <- ur[rownames(end.cor.integrated@assays$integrated@data),]
ar <- end.cor.integrated@assays$RNA@counts
ar <- ar[rownames(end.cor.integrated@assays$integrated@data),]
sr <- sr/ar;
ur <- ur/ar;
sr <- as.matrix(sr)
ur <- as.matrix(ur)
sr[is.nan(sr)] = 0;
ur[is.nan(ur)] = 0;
colnames(sr) <- colnames(end.cor.integrated)
colnames(ur) <- colnames(end.cor.integrated)
int <- as.matrix(end.cor.integrated@assays$integrated@counts)
#int[which(int < 0)]=0
spliced <- int*sr;
unspliced <- int*ur;
sg <- intersect(rownames(spliced), rownames(unspliced));
spliced <- spliced[match(sg, rownames(spliced)),];
unspliced <- unspliced[match(sg, rownames(unspliced)),];
meta <- end.cor.integrated@meta.data[,grep("time.anno|celltype.anno|time.celltype.anno|CytoTRACE",colnames(end.cor.integrated@meta.data))];
var <- sg;
pca_int <- end.cor.integrated@reductions$pca@cell.embeddings;
umap_int <- end.cor.integrated@reductions$umap@cell.embeddings;
save(spliced, unspliced, meta, var, pca_int, umap_int, file = "./supp_data/Epidermis_LRC_Atlas_scVelo_input.RData")
end.cor.traj.idx <- which(rc.integrated$celltype.anno == "Quiescent Center" | rc.integrated$celltype.anno == "Columella");
end.cor.integrated <- subset(rc.integrated, cells = colnames(rc.integrated)[end.cor.traj.idx])
end.cor.integrated <- RunUMAP(end.cor.integrated, reduction = "pca", dims = 1:8, umap.method = "umap-learn", metric = "correlation")
plot_anno_c <- function(rc.integrated){
order <- c("Quiescent Center","Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Metaphloem & Companion Cell","Protophloem", "Xylem", "Procambium","Xylem Pole Pericycle","Phloem Pole Pericycle", "Protoxylem", "Metaxylem", "Unknown")
palette <- c("#9400d3","#5ab953", "#bfef45", "#008080", "#21B6A8", "#82b6ff", "#0000FF","#e6194b", "#dd77ec", "#9a6324", "#ffe119", "#ff9900", "#ffd4e3", "#9a6324", "#ddaa6f", "#EEEEEE")
rc.integrated$celltype.anno <- factor(rc.integrated$celltype.anno, levels = order[sort(match(unique(rc.integrated$celltype.anno),order))])
color <- palette[sort(match(unique(rc.integrated$celltype.anno),order))]
p1 <- DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno", cols=color)
p2 <- DimPlot(rc.integrated, reduction = "umap", group.by = "time.anno",order = c("Distal Columella","Proximal Columella","Meristem"),cols = c("#e5f5e0", '#fee0d2','#de2d26'))
options(repr.plot.width=20, repr.plot.height=8)
gl <- lapply(list(p1, p2), ggplotGrob)
gwidth <- do.call(unit.pmax, lapply(gl, "[[", "widths"))
gl <- lapply(gl, "[[<-", "widths", value = gwidth)
gridExtra::grid.arrange(grobs=gl, ncol=2)
}
end.cor.integrated[["umap"]]@cell.embeddings[,1] <- end.cor.integrated[["umap"]]@cell.embeddings[,1]*-1
end.cor.integrated[["umap"]]@cell.embeddings[,2] <- end.cor.integrated[["umap"]]@cell.embeddings[,2]*-1
plot_anno_c(end.cor.integrated)
expression_matrix <- end.cor.integrated@assays$integrated@data
expression_matrix[which(expression_matrix < 0)]=0
expression_matrix <- as(expression_matrix, "dgCMatrix")
end.cor.integrated@assays$integrated@counts <- expression_matrix
results <- CytoTRACE(as.matrix(end.cor.integrated@assays$integrated@counts), ncores = 16, subsamplesize = 1000)
end.cor.integrated$CytoTRACE <- results$CytoTRACE
The number of cells in your dataset exceeds 3,000. CytoTRACE will now be run in fast mode (see documentation). You can multi-thread this run using the 'ncores' flag. To disable fast mode, please indicate 'enableFast = FALSE'. Warning message in CytoTRACE(as.matrix(end.cor.integrated@assays$integrated@counts), : “8 genes have zero expression in the matrix and were filtered” CytoTRACE will be run on 9 sub-sample(s) of approximately 966 cells each using 9 / 16 core(s) Pre-processing data and generating similarity matrix... Warning message in CytoTRACE(as.matrix(end.cor.integrated@assays$integrated@counts), : “1 poor quality cells were filtered based on low or no expression. See 'filteredCells' in returned object for names of filtered cells.” Calculating gene counts signature... Smoothing values with NNLS regression and diffusion... Calculating genes associated with CytoTRACE... Done
options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(end.cor.integrated, features = "CytoTRACE", pt.size=0.5)+ scale_colour_gradientn(colours = rev(brewer.pal(11,"Spectral")))
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
saveRDS(end.cor.integrated,'./supp_data/Columella_Atlas.rds')
# Prepare scVelo input
rc.su.counts <- readRDS('./Root_Atlas_spliced_unspliced_raw_counts.rds')
rc.su.counts <- subset(rc.su.counts, cells = colnames(rc.su.counts)[end.cor.traj.idx])
sr <- rc.su.counts@assays$spliced_RNA@counts
sr <- sr[rownames(end.cor.integrated@assays$integrated@data),]
ur <- rc.su.counts@assays$unspliced_RNA@counts
ur <- ur[rownames(end.cor.integrated@assays$integrated@data),]
ar <- end.cor.integrated@assays$RNA@counts
ar <- ar[rownames(end.cor.integrated@assays$integrated@data),]
sr <- sr/ar;
ur <- ur/ar;
sr <- as.matrix(sr)
ur <- as.matrix(ur)
sr[is.nan(sr)] = 0;
ur[is.nan(ur)] = 0;
colnames(sr) <- colnames(end.cor.integrated)
colnames(ur) <- colnames(end.cor.integrated)
int <- as.matrix(end.cor.integrated@assays$integrated@counts)
#int[which(int < 0)]=0
spliced <- int*sr;
unspliced <- int*ur;
sg <- intersect(rownames(spliced), rownames(unspliced));
spliced <- spliced[match(sg, rownames(spliced)),];
unspliced <- unspliced[match(sg, rownames(unspliced)),];
meta <- end.cor.integrated@meta.data[,grep("time.anno|celltype.anno|time.celltype.anno|CytoTRACE",colnames(end.cor.integrated@meta.data))];
var <- sg;
pca_int <- end.cor.integrated@reductions$pca@cell.embeddings;
umap_int <- end.cor.integrated@reductions$umap@cell.embeddings;
save(spliced, unspliced, meta, var, pca_int, umap_int, file = "./supp_data/Columella_Atlas_scVelo_input.RData")
end.cor.traj.idx <- which(rc.integrated$celltype.anno == "Xylem Pole Pericycle" | rc.integrated$celltype.anno == "Phloem Pole Pericycle" | rc.integrated$celltype.anno == "Metaxylem" | rc.integrated$celltype.anno == "Protoxylem" | rc.integrated$celltype.anno == "Procambium" | rc.integrated$celltype.anno == "Quiescent Center"| rc.integrated$celltype.anno == "Metaphloem & Companion Cell"| rc.integrated$celltype.anno == "Protophloem");
end.cor.integrated <- subset(rc.integrated, cells = colnames(rc.integrated)[end.cor.traj.idx])
end.cor.integrated <- RunUMAP(end.cor.integrated, reduction = "pca", dims = 1:50, umap.method = "umap-learn", metric = "correlation")
end.cor.integrated[["umap"]]@cell.embeddings[,1] <- end.cor.integrated[["umap"]]@cell.embeddings[,1]*-1
u2 <- end.cor.integrated@reductions$umap@cell.embeddings[,1]
u1 <- end.cor.integrated@reductions$umap@cell.embeddings[,2]
end.cor.integrated@reductions$umap@cell.embeddings[,1] <- u1
end.cor.integrated@reductions$umap@cell.embeddings[,2] <- u2
plot_anno(end.cor.integrated)
expression_matrix <- end.cor.integrated@assays$integrated@data
expression_matrix[which(expression_matrix < 0)]=0
expression_matrix <- as(expression_matrix, "dgCMatrix")
end.cor.integrated@assays$integrated@counts <- expression_matrix
results <- CytoTRACE(as.matrix(end.cor.integrated@assays$integrated@counts), ncores = 16, subsamplesize = 1000)
end.cor.integrated$CytoTRACE <- results$CytoTRACE
The number of cells in your dataset exceeds 3,000. CytoTRACE will now be run in fast mode (see documentation). You can multi-thread this run using the 'ncores' flag. To disable fast mode, please indicate 'enableFast = FALSE'. CytoTRACE will be run on 35 sub-sample(s) of approximately 1009 cells each using 16 / 16 core(s) Pre-processing data and generating similarity matrix... Calculating gene counts signature... Smoothing values with NNLS regression and diffusion... Calculating genes associated with CytoTRACE... Done
options(repr.plot.width=8, repr.plot.height=8)
FeaturePlot(end.cor.integrated, features = "CytoTRACE", pt.size=0.5)+ scale_colour_gradientn(colours = rev(brewer.pal(11,"Spectral")))
Scale for 'colour' is already present. Adding another scale for 'colour', which will replace the existing scale.
saveRDS(end.cor.integrated,'./supp_data/Stele_Atlas.rds')
# Prepare scVelo input
rc.su.counts <- readRDS('./Root_Atlas_spliced_unspliced_raw_counts.rds')
rc.su.counts <- subset(rc.su.counts, cells = colnames(rc.su.counts)[end.cor.traj.idx])
sr <- rc.su.counts@assays$spliced_RNA@counts
sr <- sr[rownames(end.cor.integrated@assays$integrated@data),]
ur <- rc.su.counts@assays$unspliced_RNA@counts
ur <- ur[rownames(end.cor.integrated@assays$integrated@data),]
ar <- end.cor.integrated@assays$RNA@counts
ar <- ar[rownames(end.cor.integrated@assays$integrated@data),]
sr <- sr/ar;
ur <- ur/ar;
sr <- as.matrix(sr)
ur <- as.matrix(ur)
sr[is.nan(sr)] = 0;
ur[is.nan(ur)] = 0;
colnames(sr) <- colnames(end.cor.integrated)
colnames(ur) <- colnames(end.cor.integrated)
int <- as.matrix(end.cor.integrated@assays$integrated@counts)
#int[which(int < 0)]=0
spliced <- int*sr;
unspliced <- int*ur;
sg <- intersect(rownames(spliced), rownames(unspliced));
spliced <- spliced[match(sg, rownames(spliced)),];
unspliced <- unspliced[match(sg, rownames(unspliced)),];
meta <- end.cor.integrated@meta.data[,grep("time.anno|celltype.anno|time.celltype.anno|CytoTRACE",colnames(end.cor.integrated@meta.data))];
var <- sg;
pca_int <- end.cor.integrated@reductions$pca@cell.embeddings;
umap_int <- end.cor.integrated@reductions$umap@cell.embeddings;
save(spliced, unspliced, meta, var, pca_int, umap_int, file = "./supp_data/Stele_Atlas_scVelo_input.RData")